Transcriptomic Analysis Revealed Key Defense Genes and Signaling Pathways Mediated by the Arabidopsis thaliana Gene SAD2 in Response to Infection with Pseudomonas syringae pv. Tomato DC3000

Nucleocytoplasmic transport receptors play key roles in the nuclear translocation of disease resistance proteins, but the associated mechanisms remain unclear. The Arabidopsis thaliana gene SAD2 encodes an importin β-like protein. A transgenic Arabidopsis line overexpressing SAD2 (OESAD2/Col-0) showed obvious resistance to Pseudomonas syringae pv. tomato DC3000 (Pst DC3000) compared to the wild type (Col-0), but the knockout mutant sad2-5 was susceptible. Transcriptomic analysis was then performed on Col-0, OESAD2/Col-0, and sad2-5 leaves at 0, 1, 2, and 3 days post-inoculation with Pst DC3000. A total of 1825 differentially expressed genes (DEGs) were identified as putative biotic stress defense genes regulated by SAD2, 45 of which overlapped between the SAD2 knockout and overexpression datasets. Gene Ontology (GO) analysis indicated that the DEGs were broadly involved in single-organism cellular metabolic processes and in response to stimulatory stress. Kyoto Encyclopedia of Genes and Genomes (KEGG) biochemical pathway analysis revealed that many of the DEGs were associated with the biosynthesis of flavonoids and other specialized metabolites. Transcription factor analysis showed that a large number of ERF/AP2, MYB, and bHLH transcription factors were involved in SAD2-mediated plant disease resistance. These results provide a basis for future exploration of the molecular mechanisms associated with SAD2-mediated disease resistance and establish a set of key candidate disease resistance genes.


Introduction
Plant pathogens are a serious threat to global food production, reducing crop yield and quality [1]. With increases in the human population and decreases in available cultivated land, methods of improving grain yield and quality have become a key focus of scientific research [2,3]. Previous studies have revealed plant disease resistance (R) genes and elucidated molecular mechanisms associated with plant disease resistance [4,5]. This research not only provides a basis for fundamental understanding of plant pathogen defense mechanisms but also lays a theoretical foundation for the more efficient improvement of crop yield and quality traits.
Because they are immobile, plants cannot avoid unfavorable environments or hazards during growth and development; they have therefore developed complex defense mechanisms [6,7]. To survive, plants must effectively sense and defend against invading pathogens [8]. To this end, they have evolved two primary layers of defense. The first layer is pathogen-associated molecular pattern (PAMP)-triggered immunity (PTI) [6,9]. PTI is triggered when PAMPs or microbe-associated molecular patterns (MAMPs) are sensed by pattern recognition receptors (PRRs) on the outer surface of plant cells. PRRs bind to co-receptors to activate a series of downstream immune responses, which initiate defense responses. Such responses include reactive oxygen species (ROS) bursts, mitogenactivated protein kinase (MAPK) activation, and defense-related gene expression. PTI can effectively prevent the invasion of most pathogenic bacteria, but the associated defense responses are typically short-lived [10,11]. To circumvent plant immune responses, pathogens produce effector proteins that are delivered directly into plant cells, bypassing PTI and making plants susceptible to disease [12,13]. In response, plants have acquired additional intracellular receptors and sensors called R proteins. These proteins recognize pathogen effectors to initiate a second layer of defense called effector-triggered immunity (ETI). ETI responses are generally stronger and more persistent than PTI responses and often lead to localized cell death [14,15]. Recognition of pathogen effectors by R proteins triggers nuclear translocation of disease-resistance proteins, thereby activating defense signals [16,17]. Nucleocytoplasmic transport receptors play key roles in this process.
In plants, a variety of nucleocytoplasmic transport pathways, including nuclear export of mRNA and nuclear import of immune-related proteins, are involved in innate immunity. Many macromolecules involved in defense responses, such as RNA and proteins, must be transferred between the cytoplasm and the nucleus to properly perform their functions. Nucleocytoplasmic transport mediates key processes: immune receptor activation; defense signal generation; intracellular defense protein transport to pathogen inoculation sites; and stimulation of programmed cell death [18][19][20]. An increasing number of studies therefore indicate that nucleocytoplasmic transport receptors play key roles in pathogen resistance [21][22][23][24][25]. Importin β is a nuclear import receptor that was first discovered in human cells; it forms a complex with importin α, which carries a substrate, to import the substrate into the nucleus [10]. Subsequently, multiple nuclear transport receptors that are homologous to importin β have been identified and classified as importin β-like proteins. Many functional disease-resistance proteins are transported to the nucleus with the assistance of importin β family proteins to induce defense responses, meaning that importin β family proteins play important roles in defense against disease [26][27][28]. Mutating the Arabidopsis gene encoding the importin β superfamily member MOS14 causes the splicing of two R genes, SNC1 and RPS4, thus impairing SNC1-and RPS4-mediated disease resistance. These findings indicate that the nuclear import of MOS14 is required for the correct splicing of SNC1 and RPS4 [26]. In Arabidopsis thaliana, SAD2 is encoded by the importinβ-like gene AT2G31660 and contains a canonical importin β domain [29][30][31]. In a previous study, we screened an Arabidopsis mutant, sad2-5, which harbors a transfer (T)-DNA insertion in the 18th exon of SAD2. In response to fumonisin-B1 (FB1) treatment, sad2-5 showed a clear cell death phenotype and abnormal H 2 O 2 accumulation. These results suggested that the nuclear transport protein SAD2 plays an important role in calcium and H 2 O 2 -mediated programmed cell death [32].
In the present study, we found that a transgenic Arabidopsis line overexpressing SAD2 (OESAD2/Col-0) was resistant to Pseudomonas syringae pv. tomato (Pst) DC3000 infection, whereas the SAD2 knockout mutant (sad2-5) was susceptible. This indicated that SAD2 may mediate the bacterial pathogen defense response in Arabidopsis. To identify candidate genes that are regulated by SAD2 and involved in response to Pst DC3000, we analyzed the transcriptomes of wild-type (Col-0), sad2-5, and OESAD2/Col-0 plants at multiple time points after treatment with Pst DC3000 to establish differentially expressed genes (DEGs) between the genotypes. Through analysis of enriched Gene Ontology (GO) functional annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways among the DEGs, we identified key stress-response genes regulated by SAD2, metabolic pathways associated with disease resistance, and disease-resistance signal transduction pathways. This study provides a theoretical reference for further exploration of resistance genes involved in the nuclear cytoplasmic transport pathway and lays the foundation for genetic improvement of crop plants through the introgression of the candidate genes identified here.

SAD2 Expression Increased Resistance to Pst DC3000 Infection
In a previous study, we found that the nuclear transport protein SAD2 was involved in calcium and H 2 O 2 -mediated programmed cell death in Arabidopsis [23]. To investigate whether SAD2 was also involved in plant pathogen responses, we here measured pathogen susceptibility in Col-0, sad2-5, and OESAD2/Col-0 plants at 4 days post-inoculation (dpi) with Pst DC3000. The sad2-5 mutant was more susceptible to Pst DC3000 infection than Col-0, whereas OESAD2/Col-0 was more resistant ( Figure 1A). Quantitative reverse transcription RT-qPCR showed that SAD2 was knocked out in sad2-5 and expressed 21.03-fold higher in OESAD2/Col-0 than in Col-0 plants ( Figure 1B). These results indicated that SAD2 may have a positive role in regulating the Arabidopsis pathogen response.

RNA-Sequencing Library Construction, Sequencing, and Read Mapping
To identify the genes regulated by SAD2 in the pathogen resistance response, Col-0, sad2-5, and OESAD2/Col-0 leaves were collected at 0, 1, 2, and 3 dpi with Pst DC3000 and analyzed via RNA-sequencing ( Table 1). The number of raw reads per sample ranged from 27,256,874 to 44,222,634, and the number of raw bases ranged from 4.08 to 6.63 Gb.
We then performed quality control on the raw sequencing data, which included removing low-quality reads, deleting adapters, calculating the sequencing error rate, calculating the Q20 and Q30 values, and determining the GC content. The number of clean reads per sample ranged from 26,714,782 to 42,603,821, and the number of clean bases ranged from 4.01 to 6.39 Gb. The per-sample mapping rate to the reference genome ranged from 96.55 to 98.17%, and the total mapping rate showed that more than 96.55% of the clean reads were successfully aligned to the reference genome. Each species has a typical percentage of GC content, which can be compared to the GC content of a given sample to indicate whether there was interference from other species [33]. The GC content here ranged from 44.75% to 46.15%, demonstrating that the samples were relatively stable and did not have interference from other species. Base quality is an important indicator for measuring sequencing accuracy. Q20 and Q30 represent base misrecognition probabilities of 1% and 0.1%, respectively. The Q30 values for our data ranged from 93.36 to 94.58%, indicating high sequencing accuracy. Overall, the data quality assessments revealed that the sequencing data were of high quality and therefore suitable for further bioinformatics analyses. We next analyzed the sequencing data to identify DEGs in the SAD2 knockout (sad2-5) and overexpression (OESAD2/Col-0) lines compared to Col-0 (Supplementary Table S1). At 0 dpi, there were a total of 994 DEGs in sad2-5 compared to Col-0, with 357 up-regulated and 637 down-regulated genes ( Figure 2A). There were 159 DEGs (62 up-regulated and 97 downregulated) at 1 dpi, 1761 DEGs (186 up-regulated and 575 down-regulated) at 2 dpi, and 82 DEGs (35 up-regulated and 47 down-regulated) at 3 dpi ( Figure 2A). In OESAD2/Col-0 samples, there were a total of 3067 DEGs at 0 dpi, with 1195 up-regulated and 1872 downregulated ( Figure 2B). There were 988 DEGs (352 up-regulated and 636 down-regulated) at 1 dpi, 471 DEGs (91 up-regulated and 380 down-regulated) at 2 dpi, and 150 DEGs (42 up-regulated and 108 down-regulated) at 3 dpi. The two mutant lines differed with respect to the dpi with the largest number of DEGs; in sad2-5 compared to Col-0, it was at 2 dpi, whereas it was at 1 dpi for OESAD2/Col-0. Further analyses were then conducted on the DEGs.

Key DEGs Regulated by SAD2 and Involved in the Pst DC3000 Defense Response
A pathogen resistance test showed that sad2-5 was more susceptible to Pst DC3000 infection than Col-0, whereas OESAD2/Col-0 exhibited resistance. This indicated that SAD2 positively regulated the pathogen defense response. Therefore, DEGs associated with pathogen resistance were expected to display contrasting expression profiles between sad2-5 and OESAD2/Col-0 samples during Pst DC3000 infection. To identify the key DEGs that were regulated by SAD2 and involved in response to Pst DC3000, we performed Venn and heatmap analyses ( Figure 3).
In sad2-5 plants, a total of 236 unique DEGs were up-regulated by Pst DC3000 over the three dpi ( Figure 3A); in OESAD2/Col-0 plants, 938 DEGs were down-regulated by Pst DC3000 ( Figure 3B). We then analyzed the genes that were common between these two sets, i.e., those that were up-regulated by Pst DC3000 infection in sad2-5 but down-regulated in OESAD2/Col-0 ( Figure 3C). This yielded a set of 24 genes in common, which were likely regulated by SAD2. Similarly, we analyzed the unique DEGs that were down-regulated by Pst DC3000 infection in sad2-5 plants (661 DEGs) ( Figure 3D), those that were up-regulated by Pst DC3000 infection in OESAD2/Col-0 plants (431 DEGs) ( Figure 3E), and the overlap between those two sets (23 genes) ( Figure 3F).

Identification of Transcription Factors (TFs) Involved in the Pathogen Response
TFs play important roles in many critical life activities in plants, including disease defense responses [34,35]. Most TFs must be transported to the nucleus to perform their functions. SAD2 belongs to a class of cytoplasmic transport receptor proteins and is responsible for cytoplasmic protein transport. We therefore focused on regulatory changes to TFs (i.e., TF expression levels) in response to SAD2 perturbation. Of the 1825 DEGs that we identified as responsive to Pst DC3000 infection, 137 DEGs were annotated as having TF activity; most of these belonged to eight main TF families: ERF/AP2, MYB, bHLH, WRKY, C2H2, GATA, NAC, and MADs-box ( Figure 5A; Supplementary Table S5). The largest number of DEGs belonged to the ERF/AP2 TF family (23), followed by the MYB family (16) and the bHLH family (10). We next analyzed the expression levels of DEGs belonging to the ERF/AP2, MYB, and bHLH families in Col-0, sad2-5, and OESAD2/Col-0 plants ( Figure 5B; Supplementary Table S6). Different expression patterns were observed for these genes between the three time points (1, 2, and 3 dpi). In total, 24 DEGs were up-regulated, and 25 were down-regulated after pathogen treatment. Notably, one set of genes that clustered together based on expression levels (AT5G07310, AT1G13300, AT5G10570, and AT2G22750) were significantly up-regulated in OESAD2/Col-0 plants at 3 dpi. Expression levels of genes in another cluster (AT4G38620, AT1G80580, AT2G40340, AT5G61890, AT3G06490, AT1G06160, AT5G61590, AT1G66390, AT4G21440, and AT4G09820) were significantly upregulated in sad2-5 plants at 3 dpi; the clustered genes AT1G26945, AT4G20970, AT2G32460, AT1G43160, and AT1G57560 were significantly up-regulated in OESAD2/Col-0 plants at 1 dpi. The results suggested that these differentially expressed TFs may have been involved in the disease defense response mediated by the nucleocytoplasmic transport receptor SAD2 in Arabidopsis.

Verification of RNA-Seq Results via RT-qPCR
To confirm the reliability of the transcriptome data, we randomly selected 10 DEGs for verification with quantitative real-time polymerase chain reaction (RT-qPCR). The gene expression patterns observed in the RT-qPCR and RNA-seq data showed similar trends (Figure 6), indicating that the transcriptome sequencing results were reliable.

Discussion
We here studied Col-0 plants, the SAD2 knockout mutant sad2-5, and the SAD2overexpression lines OESAD2/Col-0 to identify defense response genes regulated by the nucleocytoplasmic transport receptor SAD2. At 4 dpi with Pst DC3000, sad2-5 leaves were more susceptible than Col-0 leaves were, whereas OESAD2/Col-0 plants exhibited a distinct disease-resistant phenotype. Through transcriptomic analysis of inoculated leaves at 0, 1, 2, and 3 dpi, we identified 45 significant DEGs shared between sad2-5 and OESAD2/Col-0 compared to Col-0 plants. These genes had contrasting expression patterns between sad2-5 and OESAD2/Col-0 and were therefore classified as key candidate SAD2-mediated defense-related genes. For example, AT2G38530, also known as a lipid transfer protein 2 (LTP2), has been reported to play an important role in pathogen defense. Exogenous expression of barley LTP2 in tobacco or Arabidopsis increases the tolerance of these plants to bacterial pathogens [36]. In the present study, the gene encoding LTP2 was significantly up-regulated in OESAD2/Col-0 plants (peaking at 1 dpi), whereas the expression level did not change in the susceptible genotype sad2-5. AT4G28940, a phosphorylase superfamily protein gene, was expressed at extremely low levels in sad2-5 and OESAD2/Col-0 before Pst DC3000 inoculation; after Pst DC3000 inoculation, the gene was continuously up-regulated, and expression levels in OESAD2/Col-0 exceeded levels in sad2-5 by~1.5-fold. Previous studies have shown that F-type ATPase and α-1, 4 glucan phosphorylase may be important mediator proteins in the defense response to Vibrio asiatica [37]. The Phytophthora soya avirulent effector protein Avr3b is a secreted NADH and ADP-ribose pyrophosphorylase that regulates plant immunity [38]. Expression of a GDP-L-galactose phosphorylase-like gene in Chinese wild grape species induces responses to necrotrophic pathogens and defense signaling molecules [39]. We hypothesize that AT4G28940 (the phosphorylase superfamily protein) may be involved in anti-disease defense responses regulated by the nucleocytoplasmic transport receptor SAD2. AT2G04040 is a MATE efflux family protein. This gene was significantly up-regulated in OESAD2/Col-0, with the highest expression level at 3 dpi. MATE proteins have important roles in plant disease defense responses through both positive and negative regulatory mechanisms [40,41]. For example, the Arabidopsis gene EDS5 encodes a MATE transporter homolog; an eds5 mutant is highly susceptible to infection with several different pathogens, whereas EDS5 overexpression in Arabidopsis increases pathogen resistance [42,43]. The MATE transporter DTX18 is a hydroxycinnamic acid amide (HCAA) coumaroylagmatine transporter with unique specificity that inhibits infection with Pseudomonas. In potato (Solanum tuberosum), DTX18 exports HCAAs to defend against the causal agent of late blight disease [44]. Furthermore, the MATE transporter ADS1 negatively regulates the expression of salicylic acid (SA) biosynthetic genes and PATHOGENESIS-RELATED GENE 1 (PR-1), increasing plant susceptibility to pathogens [45]. Likewise, overexpressing either one of two rice MATE transporter genes, OsMATE1 or OsMATE2, in Arabidopsis increases susceptibility to Pst DC3000 [46]. Therefore, we speculate that these genes are involved in regulating pathogen defense responses.
The 45 putative SAD2-mediated defense-related DEGs were primarily enriched in the GO terms "secondary metabolic process" (GO:0019748, 7) and "flavonoid biosynthesis process" (GO:0009813, 5). GO annotation analysis revealed that SAD2 was involved in the regulation of disease resistance through single-organism cellular metabolic processes and responses to stress stimuli. The subcellular localization annotations associated with these processes were primarily the extracellular domain, plasma membrane, and intracellular membrane. Lipids are an integral part of cell membranes and play important roles in activating and mediating signal transduction pathways [47]. DEGs were also enriched in catalytic activity, nucleic-acid-binding TF activity, and oxidoreductase activity. AT1G66100, which was significantly up-regulated in OESAD2/Col-0, was annotated as a plant thionin that was mainly present in the extracellular region (GO: 0005576) and was involved in biological processes such as "response to stimulus" (GO: 0050896). Thionin is hypothesized to induce pore formation in pathogenic cell membranes, allowing potassium and calcium ion leakage [48,49]. Numerous studies have shown that the secreted antifungal compound thionin ASTHI2.4 inhibits the toxicity of fungal fruiting body lectins in Fusarium graminearum [50]. In citrus, overexpression of a modified plant thionin enhances disease resistance to citrus canker and Huanglongbing (HLB) [51]. Transgenic sweet potatoes overexpressing thionin are resistant to black rot caused by Ceratocystis fimbriata in the leaves and storage roots [52]. Therefore, we posit that AT1G66100 (a plant thionin) was involved in the bacterial defense response regulated by SAD2.
KEGG analysis revealed that the most highly enriched biochemical pathways among the DEGs were "metabolic pathways" and "biosynthesis of secondary metabolites". In addition, the 45 putative key DEGs were significantly enriched in "biosynthesis of secondary metabolites" and "flavonoid biosynthesis", consistent with the GO term analysis.
Flavonoids are a class of plant-specialized metabolites (PSMs) that protect against pathogen infection [53]. PSMs have important functions in mediating interactions between plants and other organisms. Metabolites produced by the phenylpropane metabolic pathway (e.g., flavonoids, phenols, and lignins) are also involved in plant pathogen defense responses [54]. Flavanone 3-hydroxylase (F3H) is a key player in the flavonoid biosynthetic pathway, which can coordinate plant defense mechanisms and enhance plant tolerance to biotic and abiotic stresses. A previous study showed that excessive flavonoid accumulation can enhance ROS scavenging in Arabidopsis, improving plant antioxidant activity and drought resistance [55]. Rice lines overexpressing F3H have enhanced tolerance to bacterial leaf blight (BLB) compared to wild-type plants due to excessive accumulation of antioxidant flavonoids [56]. In addition, several members of the flavonoid biosynthetic pathway were differentially expressed in sad2-5 and OESAD2/Col-0 compared to Col-0, namely AT4G22880 (leucoanthocyanidin dioxygenase [LDOX]), AT5G13930 (TT4, a chalcone and stilbene synthase family protein), and AT5G07990 (TT7, a cytochrome P450 superfamily protein). This pathway is also closely related to plant disease resistance responses. SA, jasmonate (JA), and ethylene (ET) are important components of plant defense signaling [57]. Pathogen-triggered biosynthesis of these defense hormones and their subsequent signaling activity induces activation of a complex set of defense-related genes, leading to multiple defense responses. The DEGs identified here as related to disease resistance, and defense responses were significantly enriched in MAPK signaling [58], starch and sucrose metabolism, and fatty acid elongation pathways. These results suggest that the flavonoid and specialized metabolite pathways have important roles in SAD2-mediated disease resistance.
TFs are important factors in biotic stress responses, such as those induced by insect pests or pathogens. Numerous studies have reported the involvement of many TF families in plant pathogen defense responses. For example, the ET response factor ERF11 enhances plant resistance to bacterial pathogens [59]. The ERF TF GmERF113 increases soybean resistance to Phytophthora soja [60]. The tomato ERF TF SlERF84 enhances tolerance to drought and salt stress but negatively regulates defense responses to Pst DC3000 [61]. Overexpression of the apple (Malus domestica) gene MdERF100 in Arabidopsis enhances resistance to powdery mildew by regulating the JA and SA signaling pathways [62]. R2R3-MYB TFs can improve brown planthopper resistance in rice by regulating the phenylalanine ammonia lyase pathway [63]. In wheat (Triticum aestivum), the TF TaMYB29 is involved in the prevention of wheat stripe rust and is significantly induced by SA, abscisic acid (ABA), JA, ET, and Pst [64]. In rice (Oryza sativa), a bHLH transcriptional activator, OsbHLH6, has been identified as an enhancer of disease resistance by dynamically regulating SA and JA signaling through nuclear-cytoplasmic trafficking [65]. MdbHLH92 is involved in plant defense against powdery mildew [62]. In addition, knocking down bHLH132 in tomato stunts plant development and enhances susceptibility to Xanthomonas euvesicatoria [66]. We here identified a total of 137 DEGs with TF activity that may have been involved in SAD2-mediated defense responses; most such TFs belonged to the EFR, MYB, and bHLH families (23, 10, and 9 DEGs, respectively). Two of these genes (AT1G26945 and AT1G10586) were bHLH DNA-binding superfamily proteins and were specifically up-regulated in OESAD2/Col-0 plants at 1 and 3 dpi. We hypothesize that they may participate in defense responses that are regulated by the nuclear cytoplasmic transport receptor SAD2.
In summary, we here analyzed RNA-sequencing data from wild-type (Col-0), SAD2 knockout mutant (sad2-5), and SAD2-overexpression (OESAD2/Col-0) plants. Each genotype was inoculated with Pst DC3000 and analyzed at multiple time points. The resulting gene expression profiles revealed key candidate genes and metabolic pathways that appeared to be involved in defense responses that are regulated by the nucleocytoplasmic transport receptor SAD2. These results improve our understanding of the mechanisms by which nuclear cytoplasmic transport receptors mediate plant defense responses. Furthermore, this study provides a foundation for future exploration of the molecular mechanisms associated with SAD2-mediated disease resistance and establishes a set of strong candidate plant disease-resistance genes

Plant Materials and Growth Conditions
Arabidopsis ecotype Col-0, the SAD2 knockout mutant sad2-5, and the SAD2overexpression line OESAD2/Col-0 were used in this study. sad2-5 is a previously described T-DNA insertion line [29,30]. OESAD2/Col-0 was generated for the present study by inserting the coding region of SAD2 into the pCAMBIA1307 vector under the control of the cauliflower mosaic virus 35S promoter. Col-0 was transformed with the recombinant plasmid by Agrobacterium-mediated transformation via floral dip, and T4 homozygous lines were used for this study.
For plants used in susceptibility analyses and RNA-seq assays, seeds were sterilized and sown on 1/2× Murashige and Skoog (MS) medium containing 2.5% (w/v) sucrose and 1% (w/v) agar. At approximately 1 week of age, seedlings were transplanted into soil and raised in a growth chamber at 22 • C with a 16/8 h light/dark cycle at 60-70% relative humidity.

Arabidopsis Inoculation with Pst DC3000 and Leaf Sample Collection
Pst DC3000 was cultured on King's B (KB) solid medium with 25 µg/mL rifampicin at 28 • C for 2 d. Single colonies were cultured in KB liquid medium overnight at 28 • C with shaking, reaching an OD600 ≥ 1. The OD600 was then adjusted to 0.002 with sterile ddH 2 O. Arabidopsis rosette leaves were infiltrated with the bacterial suspension using 1 mL needleless syringes. Leaves were collected in biological triplicate at 0, 1, 2, and 3 dpi, immediately frozen in liquid nitrogen, and stored at −80 • C prior to RNA extraction.

RNA Extraction, cDNA Library Construction, and Illumina Deep Sequencing
Total RNA was extracted from leaf samples using TRIzol reagent (Invitrogen, Thermo Fisher Scientific, Beijing, China) following the manufacturer's protocol. The purity, concentration, and quality of RNA samples were measured using a NanoDrop 2000 spectrophotometer (NanoDrop, Wilmington, DE, USA) to ensure sample suitability. RNA quality was also verified through visualization of~500 ng RNA per sample on a 1% agarose gel. From suitably high-quality RNA samples, library construction and Illumina sequencing were completed at Allwegene Technology (Beijing, China). The HiSeq/NovaSeq PE150 platform was used and resulted in~6 Gb of clean data per sample.

RNA-Sequencing Data Analysis
Raw reads obtained from the sequencing facility included adaptor sequences and low-quality sequences. To ensure the accuracy of downstream analyses, we used fastp [67] to clean raw reads and FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/ fastqc, accessed on 20 October 2020) to evaluate the quality. Clean reads were aligned to the TAIR10 Arabidopsis reference genome with Hisat2 [68]. After alignment, transcript abundance was quantified for each gene with SAMtools (http://samtools.sourceforge.net/, accessed on 25 December 2020) [69] and HTSeq (http://www-huber.embl.de/HTSeq/ doc/overview.html, accessed on 8 March 2021) [70]. The Bioconductor 'DESeq2' [71] package in R was used to identify DEGs between genotypes at each time point. For gene expression analysis, the number of unambiguous tags for each gene was calculated and normalized to the number of transcripts per million tags (TPM). Genes with FDR ≤ 0.05 and |log2(FC)| ≥ 1 were classified as significant DEGs.

Gene Enrichment Analysis
TBtools [72] was used to generate Venn diagrams to display unique and overlapping DEGs between groups. The 'heatmap' package in R was used to draw heatmaps. To determine the biological functions of DEGs, we used the online tool AgriGO (http:// systemsbiology.cau.edu.cn/agriGOv2/index.php, accessed on 10 October 2021) to perform GO enrichment analysis [73]. GO terms with FDR ≤ 0.05 were considered significantly enriched. For biochemical pathway analysis, we used DAVID (https://david.ncifcrf.gov/ home.jsp, accessed on 16 November 2021) to convert each gene ID to the corresponding Entrez gene ID, then entered those gene IDs into the KOBAS website (http://kobas.cbi.pku. edu.cn/kobas3/genelist/, accessed on 16 November 2021). KEGG biochemical pathway enrichment analysis was conducted for DEGs using the KEGG database [74], with p ≤ 0.05 used as the threshold for classifying significantly enriched pathways.

Validation of RNA-Seq Data
The primers used for RT-qPCR (Supplementary Table S7) were designed with Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/, accessed on 23 October 2022) and synthesized by Sangon Biotech (Shanghai, China). ACTIN1 (AT2G37620) was used as the internal reference gene to normalize expression levels using the 2 −∆∆Ct method [75]. Approximately 2 µg of total RNA was used for first-strand cDNA synthesis (Tiangen Biotech, Beijing, China). RT-qPCR was performed in 20 µL reaction volumes using SYBR Premix Ex TaqTM II (TAKARA BIO INC